%!TeX encoding = utf-8
%!TeX program = xelatex
\documentclass[12pt,a4paper]{article}                                                                            
\usepackage{amsmath,amsfonts,amsthm,amssymb}%AMS数学符号字体等的宏包
\usepackage{amsrefs}                           %ams参考文献宏包, 不需要的用户请删除此行!
%\usepackage{hyperref}                          %生成超链接
%\usepackage[colorlinks,linkcolor=black,anchorcolor=black,citecolor=black]{hyperref}

\usepackage{geometry}                        %页面设置
\usepackage{xeCJK}                               %中文字体宏包
\usepackage{graphicx}                          %插图
\usepackage{setspace}  
\usepackage{caption}
\captionsetup{font={small}}
\usepackage{listings} %C, C++等代码高亮，需要配合xcolor宏包支持
\usepackage{tabularx}
\usepackage{xcolor}
\usepackage{fontspec}
\newfontfamily\Monaco{Monaco}

%%%%%%%%%%实用功能%%%%%%%%%%%%%%%%%%%%%%
\usepackage{indentfirst} % 解决 CJK 手段缩进问题,\setlength{\parindent}{2em}

%%%%%%%%%算法伪代码%%%%%%%%%%
\usepackage{algorithm}
\usepackage{algorithmicx}
\usepackage{algpseudocode}
\floatname{algorithm}{算法}
\renewcommand{\algorithmicrequire}{\textbf{输入:}}
\renewcommand{\algorithmicensure}{\textbf{输出:}}


\usepackage{cite}

%%%%%%%%%页面设置%%%%%%%%%%%
\setlength{\textwidth}{162mm}
\setlength{\textheight}{231mm}
\setlength{\headheight}{0cm}
\setlength{\topmargin}{-0.1cm}
\setlength{\oddsidemargin}{0cm}
\setlength{\evensidemargin}{0cm}
\setlength{\parskip}{1mm}
\setlength{\unitlength}{1mm}
\setlength{\parindent}{2.06em}
\linespread{1.5}

%%%%%%%%%%重命名%%%%%%%%%%
\renewcommand\refname{参考文献}
\renewcommand\contentsname{目录}
\renewcommand{\figurename}{图}
\renewcommand{\tablename}{表}
\renewcommand\abstractname{\bf\large 摘~~要}


%%%%%%%%%%定义编号规则%%%%%%%%%%
\renewcommand{\thefigure}{\arabic{section}.\arabic{figure}}
\renewcommand{\thetable}{\arabic{section}.\arabic{table}}
\renewcommand{\theequation}{\arabic{section}.\arabic{equation}}



%\setmonofont[Mapping={}]{Monaco}	%英文引号之类的正常显示，相当于设置英文字体
%\setsansfont{Monaco} %设置英文字体 Monaco, Consolas,  Fantasque Sans Mono
%\setmainfont{Monaco} %设置英文字体
 %\setCJKmainfont{方正兰亭黑简体}  %中文字体设置
 %\setCJKsansfont{华康少女字体} %设置中文字体
% \setCJKmonofont{华康少女字体} %设置中文字体
%%%%%%%%%%CJK字体设置%%%%%%%%%%
\setCJKmainfont[BoldFont=PingFangSC-Semibold, ItalicFont=PingFangSC-Thin]{PingFangSC-Regular}
\setCJKsansfont[BoldFont=PingFangSC-Semibold]{PingFangSC-Regular}
\setCJKmonofont{PingFangSC-Regular}

\newcommand\dd{\mathrm{d}}

%%%%%%%%%修改行间距%%%%%%%%%%%%
\usepackage{setspace}%使用间距宏包

\begin{document}

\definecolor{mygreen}{rgb}{0,0.6,0}
\definecolor{mygray}{rgb}{0.5,0.5,0.5}
\definecolor{mymauve}{rgb}{0.58,0,0.82}
% Listing环境全局设置
\lstset{
%numbers=left,
%numberstyle=\tiny,
%basicstyle=\small,
keywordstyle=\color{blue!70}, 
commentstyle=\color{red!50!green!50!blue!50},
%frame=shadowbox,
rulesepcolor=\color{red!20!green!20!blue!20},
escapeinside=``,
%xleftmargin=2em,
%xrightmargin=2em,
%aboveskip=1em
backgroundcolor=\color{white},   % choose the background color
basicstyle=\footnotesize\ttfamily,        % size of fonts used for the code
columns=fullflexible,
breaklines=true,                 % automatic line breaking only at whitespace
captionpos=b,                    % sets the caption-position to bottom
tabsize=4,
commentstyle=\color{mygreen},    % comment style
%escapeinside={\%*}{*)},          % if you want to add LaTeX within your code
keywordstyle=\color{blue},       % keyword style
stringstyle=\color{mymauve}\ttfamily,     % string literal style
frame=single,
}





\def\chntoday{\the\year~年~\the\month~月~\the\day~日}

\title{软件包HYPRE的概述与试用}
\author{姚伟鹏\footnote{电话：13370130802，邮箱：ywp@pku.edu.cn，学号：1501110118 }}
\date{\chntoday}

\maketitle                            % 生成标题

\tableofcontents                      % 插入目录

\newpage

\section{HYPRE的概述}

\subsection{HYPRE的简介}
HYPRE (High Performance Preconditioners) ，又称“高性能预条件子”，是由美国加州大学和劳伦斯—利弗莫尔国家实验室应用科学计算中心开发的，
最初用于模拟国防、环境、能源和生物科学中的现象。
如今主要用于大规模并行计算机上求解大型稀疏线性方程组——为用户提供高级并行预条件子。

下面，我们将分别介绍HYPRE的主要特征，网格接口，迭代法和预条件子以及其构造步骤等方面。

\subsection{HYPRE的主要特征}
HYPRE的主要特征为：

\begin{itemize}
\item \textbf{可扩展的预条件子}：HYPRE包含可扩展求解超大规模稀疏线性方程组的几类预条件算法，其中包括SMG (Semicoarsening MultiGrid)，又称“结构化多重网格”算法，以及基于矩阵元的AMG (Algebraic Multigrid)，又称“代数多重网格”方法。而且，在HYPRE中的BoomerAMG不但可以当做预条件子，还能直接当做求解器。
\item \textbf{常用的迭代法实现}：HYPRE提供一些最常用的基于Krylov子空间迭代法。比如求解非对称矩阵的GMRES (Generalized Minimal Residual Algorithm) 法，又称“广义最小残差算法”，以及求解对称矩阵的CG (Conjugate Gradient) 法，又称“共轭梯度”法。
\item \textbf{以网格为中心的界面}: HYPRE通过各种网格界面表示和处理稀疏矩阵，每个界面提供对一些求解器的访问，因此不需要用户去学习和创建复杂的数据结构。其中，包括基于模板的结构网格界面（Structured-Grid System interface）和半结构网格界面（Semi-Structured-Grid System interface），基于有限元的无结构界面（Unstructured interfaces）以及基于线性代数的界面。
\end{itemize}

\subsection{HYPRE的网格接口}
HYPRE库目前支持四种网格接口，它们分别为：
\begin{itemize}
\item \textbf{结构化网格接口（Struct）}：面向结构网格离散的应用。每个网格点的离散格式有相同的模式，如有限差分格式。
\item \textbf{半结构化网格接口（SStruct）}：面向半结构网格离散的应用，如局部加密 AMG，块结构网格上的应用。具体而言，包括有限差分方法和有限体积方法。
\item \textbf{有限元接口（FEI）}：面向通过有限元离散得到的线性方程。
\item \textbf{线性代数系统接口（IJ）}：传统线性代数的非结构矩阵界面。 该接口以矩阵方式显式的表示线性代数方程组，是适用范围最广泛的接口。
\end{itemize}

在HYPRE中，数据结构、求解器与网格接口的关系可以用下图表示：
\begin{figure}[!htp]
	\centering
	\includegraphics[width=0.9\textwidth]{files/interface.png}
	\label{fig1}
	\caption{其中，第一层表示的是各种线性系统的网格界面（最左边为Struct，最右边则为IJ），第二层则为各种线性求解器（或者迭代法与预条件子），第三层则表示各种数据划分和矩阵向量储存策略。}
\end{figure}

\subsection{迭代法和预条件子}

HYPRE的迭代方法如下：
\begin{itemize}
\item Krylov 解法器，如：共轭梯度法（Conjugate Gradient, CG），广义最小残差法（Generalized Minimal Residual Algorithm, GMRES）等。
\item 并行代数多重网格解法器（BoomerAMG）。
\item SuperLU直接解法器等。
\end{itemize}

HYPRE的预条件子如下：
\begin{itemize}
\item Diagonal：对角，块对角（Jacobi）预条件子。
\item PILUT：具有阙值（threshold）的并行不完全LU分解（PILU）。
\item Euclid：并行ILU预条件子的扩展。
\item SMG：并行半粗化多重网格预条件子。
\item PFMG：使用简单点光滑的并行半粗化多重网格预条件子。
\item BoomerAMG：并行代数多重网格预条件子。
\item ParaSails：并行稀疏近似逆预条件子。
\end{itemize}

在使用中，预条件子的选择和构造放在迭代法的选择和构造之前，然后再将其传递给迭代法，作为迭代法构造过程的一部分。

HYPRE为不同的接口定义了不同的数据结构，并配以适合该接口的求解器，如下图所示：
\begin{figure}[!htp]
	\centering
	\includegraphics[width=0.45\textwidth]{files/solvers.png}
	\label{fig2}
	\caption{X表示HYPRE中的某个网格界面支持某个求解器。}
\end{figure}

\subsection{求解器构造步骤}
假设在一个线性系统$Ax = b$中（$A$为系数矩阵，$b$是右端项，$x$是其解），构造一个求解器需要进行下述步骤：
\begin{lstlisting}[language={[ANSI]C},
			basicstyle=\small\Monaco]
/* 创造一个求解器，其名字假设叫SOLVER */
int HYPRE_SOLVERCreate(MPI_COMM_WORLD, &solver);
/* 设置所需要的参数 */
HYPRE_SOLVERSet(); 
// 将预条件子传递给求解器
HYPRE_SetPrecond(); 
// 这个参数用来控制误差
HYPRE_SOLVERSetTol(solver, 1.e-8);  
...
/* 构造求解器 */
// A, b, x需要提前用类似的步骤构造好
HYPRE_SOLVERSetup(solver, A, b, x); 
/* 求解线性系统 */
HYPRE_SOLVERSolve(solver, A, b, x);
/* 释放求解器 */
HYPRE_SOLVERDestroy(solver);
\end{lstlisting}


\section{稀疏矩阵}
在试用HYPRE求解具体问题之前，有必要对其最擅长的求解对象——稀疏矩阵——做一些说明及概述。

一般而言，当一个矩阵 $A\in R^{n\times n}$ 的绝大多数元素为零时，我们称其为一个稀疏矩阵。
具体的说，如果一个矩阵的非零元素的总数为 $O(n)$，则可称之为稀疏矩阵。
在微分方程数值解法中出现的线性代数方程组中，系数矩阵往往很大，但其中非零元素所占的比例却一般很小，所以稀疏矩阵是我们经常遇到的。

\subsection{稀疏矩阵的储存格式}
为了节省内存，在储存稀疏矩阵时，我们很自然的会期望只储存其中的非零元素。
然而，由于稀疏矩阵中的非零元素分布没有任何规律，所以在进行压缩存储的时侯除了需要存储非零元素的值之外，还要存储非零元素在矩阵中的位置——非零元素所在的行号和列号。

稀疏矩阵的储存格式有很多（如坐标格式（Coordinate Format），压缩稀疏行格式（Compress Sparse Row (CSR) Format），压缩稀疏列格式（Compress Sparse Col (CSC) Format），对角元优先格式（Modified Sparse Row (MSR) Format））等。

在HYPRE中使用的储存格式为CSR，对于一个稀疏矩阵A，它需要储存以下三个一维数组：
\begin{itemize}
\item \textbf{AA(i), i = 1,n}： 实型数组，按行存放矩阵 A 中的所有非零元素。
\item \textbf{JA(i), i = 1,n}： 整型数组，依次存放数组 AA 中对应元素的列号。
\item \textbf{IA(i), i = 1,n+1}： 整型数组，依次存放各行在数组 AA 或 JA 中存储的首地址。利用数组 IA，可以很方便的遍历矩阵A各行的首元素，如当 j 从 IA(i) 到 IA(i+1)-1 变化时，AA(j) 即表示矩阵 A 第 i 行的所有非零元素。
\end{itemize}


\subsection{共轭梯度法（CG）}

共轭梯度法（CG）是求解正定对称矩阵的最重要的迭代法之一，适用于对称的稀疏矩阵线性方程组。
其优点是所需存储量小，稳定性高，而且不需要任何外加参数。

用共枙梯度法求解对称正定线性方程组第 k 次迭代实质上是求 $x_k\in x_0+\mathcal{K}(A,r_0,k)$ 使得
\begin{equation*}
  \varphi(x_k) = min\{\varphi(x):x\in x_0+\mathcal{K}(A,r_0,k)\}
\end{equation*}

然而，当矩阵A并非对称正定矩阵时，$\varphi$ 不一定有极小值， 
因此不能直接用共枙梯度法求解一般线性方程组。
但这种思想却是可以推广到一般情形的。

\subsection{广义最小残差法（GMRES）}

广义最小残差法（GMRES）则是求解
$x_k\in x_0+\mathcal{K}(A,r_0,k)$, 使得
\begin{equation*}
  \| b-Ax_k\|_2 = min \{\| b-Ax\|_2: \quad x\in x_0+\mathcal{K}(A,r_0,k) \}
\end{equation*}
其优点是不受矩阵对称性的限制。

\subsection{多重网格法（MG）}

多重网格也是一种重要的迭代算法，它包含三个要素：光滑算子，限制算子和延拓算子。
一个简单的例子是：用分片线性插值作为延拓，相邻点的加权平均作为限制，松弛迭代（Guass-Seidel, Jacobi, SSOR）等简单迭代作为光滑。

多重网格解法器是HYPRE的重要特色。
它提供了多种多重网格解法器，如AMS，SMG，PFMG，BoomerAMG。
这些可以满足各种应用的需求。
其中 SMG 和 BoomerAMG 是目前实际应用中使用最广泛的两个解法器。

在下一章的算例当中，我们也会使用这种算法。


\section{HYPRE的试用}
在了解了程序包HYPRE的使用方法以及稀疏矩阵的储存与求解方法之后，现在我们来尝试使用HYPRE程序包求解一个简单的算例，
并通过比较不同迭代方法或预条件子的结果，来比较它们各自的优缺点。
最后，进行常规的并行效率分析并用HYPRE推荐使用的GLVis可视化工具给出数值结果。

\subsection{简单算例}
这里，我们选择HYPRE程序包中\textbf{src/example}下的\textbf{ex5.c}作为简单的算例，求解在$N \times N$网格上的二维拉普拉斯方程。
每个核上的网格为$n \times n$，空间分辨率为$h=1/(n+1)$。
采用IJ界面和五点差分离散格式，并考虑简单的边界条件。

将解$x$分解为$x_i$和$x_b$，它们分别代表内部解与边界解。
这样，我们就可以把矩阵A分解为：
\begin{equation*}
A = \left[
\begin{array}{cc}
A_{ii} & A_{ib} \\ 
A_{bi} & A_{bb} \\
\end{array}
\right]
\end{equation*}

将边界条件设定为简单狄利克雷边界，那么有$x_b = u_0$。
且我们只需要求解线性系统的右端项的内部解$b_i$：
\begin{equation*}
A_{ii} x_i + A_{ib} u_0= b_i 
\end{equation*}

为求简单，我们设$u_0 = 0$，那么我们只需要求解$A_{ii} x_i = b_i$。

针对这个问题，可以使用多种迭代方法，如AMG，CG，CG+AMG以及GMRES+AMG。
下面，我们就尝试使用不同的方法来对这个简单的算例进行计算。

\subsection{求解步骤}

\begin{lstlisting}[language={[ANSI]C},
			basicstyle=\footnotesize\Monaco]
/* 创造一个矩阵A，ilower和iupper限定了网格空间的左下坐标与右上坐标 */
HYPRE_IJMatrixCreate(MPI_COMM_WORLD, ilower, iupper, ilower, iupper, &A);
/* 选择并行CSR储存方式 */
HYPRE_IJMatrixSetObjectType(A, HYPRE_PARCSR);
/* 初始化A */
HYPRE_IJMatrixInitialize(A);
/* 为矩阵A的每一列赋值*/
HYPRE_IJMatrixSetValues(A, 1, &nnz, &i, cols, values);
/* 确认矩阵A已经“准备好了” */
HYPRE_IJMatrixAssemble(A);
/* 取出A以备使用 */
HYPRE_IJMatrixGetObject(A, (void**) &parcsr_A);

/* 创造右端项并赋值、确认、使用；将b换成x，则对矢量解进行同样的操作 */
HYPRE_IJVectorCreate(MPI_COMM_WORLD, ilower, iupper,&b);
HYPRE_IJVectorSetObjectType(b, HYPRE_PARCSR);
HYPRE_IJVectorInitialize(b);
HYPRE_IJVectorSetValues(b, local_size, rows, rhs_values);
HYPRE_IJVectorAssemble(b);
HYPRE_IJVectorGetObject(b, (void **) &par_b);

/* 通过参数solver_id选择迭代方法XXX来求解 */
if (solver_id == X)
...
HYPRE_XXXCreate(&solver);
HYPRE_XXXSet...();
HYPRE_XXXSetup(solver, parcsr_A, par_b, par_x);
HYPRE_XXXSolve(solver, parcsr_A, par_b, par_x);
...
/* 输出结果信息 */
HYPRE_XXXNumIterations(solver, &num_iterations);
HYPRE_XXXGetFinalRelativeResidualNorm(solver, &final_res_norm);

HYPRE_XXXDestroy(solver);
\end{lstlisting}

\subsection{结果比较与分析}
假设$N=1024, np=4,tol=1.e-6$，使用上述四种方法求解算例，所需的迭代次数与计算时间如下表所示：
\begin{table}[h]
  \small
  \begin{tabular}{|c|c|c|c|c|c|}
    \hline
          & \textbf{AMG} & \textbf{CG} & \textbf{CG+AMG} & \textbf{GMRES+AMG} & \textbf{CG+ParaSail}\\
    \hline
    迭代步数  & 9 & 1672 & 6  & 6 & 689\\ 
    \hline
    迭代时间(s) & 1.076485e+00 & 1.338007e+01 & 1.108234e+00 & 1.504573e+00 & 1.876050e+01\\      
    \hline
    迭代误差 & 3.941944e-07 & 9.953249e-07 & 7.675342e-07 & 4.949916e-08 & 9.743934e-07 \\      
    \hline
  \end{tabular}
\end{table}

从上表不难得出：
\begin{itemize}
\item 以迭代步数这一指标来看，AMG方法最好——一方面，与CG法相比，降低了100倍之多！另一方面，其他迭代步数较多的方法，如CG，在使用了AMG预条件子之后，达到一定误差所需的迭代步数得到了极大的优化；进一步，CG法在使用AMG和ParaSail两种不同的预条件子优化后，AMG的优化效果明显要高于ParaSail。
\item 虽然使用AMG 方法的迭代步数仅仅需要9步就能达到控制误差，而CG法需要1672步，但是从时间来看，使用AMG所需的时间并没有降低相应的倍数，仅仅10倍。这说明CG法的单步速度相比AMG而言，要快很多。这与它所需的储存量小，稳定性高的优点有关。 
\end{itemize}

\subsection{并行效率}
接下来，我们使用AMG迭代法，在$N=4096$的情况下，控制误差在$1.0e-6$之内，分析其并行效率。如下表所示：

\begin{table}[h]
  \begin{tabular}{|c|c|c|c|c|}
    \hline
   进程数 &1&2&3&4 \\
    \hline
    运行时间(s) & 4.077742e+01 & 3.919587e+01 & 3.418898e+01 & 3.380549e+01 \\
    \hline
  \end{tabular}
\end{table}

不难看出，此时HYPRE的并行效率并不高。可能的原因为我们尝试的算例规模过小——当规模较小时，HYPRE并不能发挥出它高效的作用。

\subsection{数值结果}

在使用AMG迭代法，$N=1024$并控制误差在$1.0e-6$之内的条件下，使用HYPRE推荐的可视化工具GLVis处理运算结果。

首先，可视化工具GLVis的安装过程比较繁琐，需要配套安装的其他两个程序包（metis与mfem），在安装时需要进行的自定义配置也较多，这可以算作GLVis的缺陷之一。
然而，相比较Matlab或者IDL而言，安装好GLVis之后，其优点无疑是大于缺陷的（求其是在处理HYPRE的数据结果的时候）。

一方面，GLVis与HYPRE很好的搭配在一起，只需要在使用HYPRE包时，加入\textbf{-vis}参数，即可由GLVis输出数据，通过HYPRE的\textbf{/src/example/vis}中提供的脚本，非常方便即可完成画图过程。
另一方面，轻量级、画图快速、对图像的操作快捷也是GLVis相对于Maltab等传统工作在作图（尤其是3D情况）方面的优点——它可以一键完成对图像的旋转、透视、切割等复杂操作，可视化效果也要比Matlab等工具所作出的三维图像好很多。

具体如下图所示：
\begin{figure}[!htp]
	\centering
	\includegraphics[width=0.3\textwidth]{files/2d.png}
%	\includegraphics[width=0.3\textwidth]{files/3d.png}
	\includegraphics[width=0.3\textwidth]{files/trans1.png} \\
	\includegraphics[width=0.3\textwidth]{files/3d_cut.png} 
%	\includegraphics[width=0.4\textwidth]{files/trans1.png}
	\includegraphics[width=0.3\textwidth]{files/trans2.png}
	\label{fig3}
	\caption{二维与三维结果图，以及GLVis对图像的切割和透明处理示意图。}
\end{figure}


\section{总结}
通过这次对HYPRE的试用作业，我学习到了很多很多。

（1）首先是对稀疏矩阵有了一个系统的了解和认识，理解了其CSR储存方式

（2）其次是认识了HYPRE这款强大的软件包，从它的下载、配置、安装，到对其手册的阅读、相关资料的查阅，最后到对其给出算例的试用，系统的学习到了如何从零开始学习使用一套软件包的方法。对其接口、使用步骤，甚至于其配套的可视化工具，都有了一个全面的认识。相信以后在再一次接触到其他软件包后，能够很快的上手使用。

（3）对于我个人来说，HYPRE中各种对线性方程求解的迭代方法和预条件子的各种参数设置是一个难点，因为它们分别是如何对求解线性方程优化的原理我还并不了解，仅仅是将其看做一个黑箱来使用，而且还会由于对其参数的设置不合理，而达不到其原本的效果。

（4）另外，由于本次算例是由HYPRE中的算例提供的，所以接口问题事实上并不存在。如果处理一个给定矩阵信息后的稀疏矩阵求解问题，那么如何根据遇到的具体问题，设置ilower和iupper、分配并行核、读取数据并向创建并初始化后的矩阵赋值，都将是潜在的问题。

（5）最后，通过算例得到的HYPRE并行效率并不高，虽然分析的可能原因是稀疏矩阵规模不够，但真正的原因是什么？还有哪些影响其效率的因素（比如对预条件子的参数设置并不合理等）？这也是一个没有弄懂的问题。


\begin{thebibliography}{99}
\bibitem{pre1} 李荣华，刘播，《微分方程数值解法(第四版)》，高等教育出版社，2009。
\bibitem{pre2} 张林波，迟学斌，莫则尧，李若，《并行计算导论》，清华大学出版社，2006。
\bibitem{pre3} 王武，迟学斌，程强，冯仰德。Hypre：高性能预条件子，超级计算通讯， 2006。
\bibitem{pre4} HYPRE Reference Manual
\bibitem{pre5} HYPRE User Manual
\end{thebibliography}









\end{document}


